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ABSTRACT 

It has been shown observationally that star formation (SF) correlates tightly with the 
presence of molecular hydrogen (H2). Therefore it would be important to investigate 
its implication on galaxy formation in a cosmological context. In the present work, we 
track the H 2 mass fraction within our cosmological smoothed particle hydrodynamics 
(SPH) code GADGET-3 using an equilibrium analytic model by Krumholz et al. This 
model allows us to regulate the star formation in our simulation by the local abundance 
of H2 rather than the total cold gas density, and naturally introduce the dependence 
of star formation on metallicity. We investigate implications of the equilibrium H 2 - 
based SF model on galaxy population properties, such as the stellar-to-halo mass 
ratio (SHMR), baryon fraction, cosmic star formation rate density (SFRD), galaxy 
specific SFR, galaxy stellar mass functions (GSMF), and Kcnnicutt-Schmidt (KS) 
relationship. The advantage of our work over the previous ones is having a large sample 
of simulated galaxies in a cosmological volume from high-redshift to z = 0. We find 
that low-mass halos with Mdm < 1O 1O - 5 M are less efficient in producing stars in the 
H 2 -based SF model at z 6, which brings the simulations to a better agreement with 
observational estimates of SHMR and GSMF at the low-mass end. This is particularly 
evident by a reduction in the number of low-mass galaxies at M* ^ 10 8 M Q in the 
GSMF. The overall SFRD is also reduced at high-z in the H 2 run, which results in 
slightly higher SFRD at low-redshift due to more abundant gas available for star 
formation at later times. This new H 2 model is able to reproduce the empirical KS 
relationship at z — naturally without the need for setting its normalization by hand, 
and overall it seems to have more advantages than the previous pressure-based SF 
model. 



1 INTRODUCTION 

Properly modeling star formation and feedback within sim- 
ulations of galaxy formation is one of the holy grails for 
computational astrophysicists. Unfortunately, current cos- 
mological simulations lack the spatial and mass resolutions 
to properly resolve the small scale processes which govern 
star formation within the interstellar medium (ISM). This 
computational restriction gives rise to the need for sub-grid 
models that can accurately describe global properties of the 
ISM. Simulation results can vary drastically depending on 
the details adopted for such sub-grid models and their feed- 
back prescriptions. It is for this reason that these sub-grid 
models rely heavily on observed empirical star formation 
models. 

The most well-known empirical star formation relation 
is the Schmidt (1959) and Kennicutt (1998) relationship, 
which relates the density (or surface density) of star for- 
mation to the gas density (or surface density), respectively. 
For numerical simulations of galaxy formation, the Schmidt 



relationship is easier to implement (e.g., Katz 1992; Cen & 
Ostriker 1992), however, observationally the Kennicutt rela- 
tionship is easier to measure because observations are done 
in the projected 2-dimensional plane. In the present work, 
we are implementing the Schmidt relationship as part of our 
SF model, but when comparing to the observations we use 
the Kennicutt relationship, hence referring to them collec- 
tively Kennicutt-Schmidt (KS) relationship. 

Recent observational evidence has suggested that star 
formation is more tightly correlated with the presence of 
molecular hydrogen (H 2 ), rather than neutral atomic (Hi) 
hydrogen (Wong & Blitz 2002; Kennicutt et al. 2007; Leroy 
et al. 2008; Bigiel et al. 2008; Bolatto et al. 2011). In partic- 
ular, Bigiel et al. (2008) studied the KS relation for a sample 
of nearby galaxies, and found little to no correlation between 
Ehi and £*, whereas £h2 was found to correlate strongly 
with Bolatto et al. (2011) used Spitzer dust continuum 
observations of the low metallicity SMC to calculate H2 sur- 
face densities without the need for a CO luminosity conver- 
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3.60 x 10 7 
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Y 
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N600L100 


100.00 


2 x 600 3 


2.78 x 10 8 


5.65 x 10 7 
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0.00 


0.00 


Y 
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Table 1. Simulation parameters used in this work. The first three simulations were used to perform tests of the H2 model and resolution 
study (Section 3.6). The second set of three simulations are the main production runs used to compare with previous SF models. The 
quantities m<j m & m gas are the particle masses of dark matter and gas particles, e is the comoving gravitational softening length, and 
z cnt j is the ending redshift of each simulation. The H2 simulations (along with N144L10 Fiducial) use an optically-thick ultra-violet 
threshold or 'OTUV (see Section 2.3; Nagamine et al. 2010). 



sion factor. Their findings suggested that H2 can be used to 
infer star formation activity even in low metallicity galaxies. 

Driven by these observational findings, new models have 
been developed relating SFRs directly to the abundance 
of H2. Some are in the form of analytic models (Fu et al. 
2010; Krumholz et al. 2008, 2009; McKee & Krumholz 2010; 
Krumholz et al. 2012), while others in the form of non- 
equilibrium, fully time-dependent calculations (Gnedin et al. 
2009; Feldmann et al. 2011; Mac Low & Glover 2012). How- 
ever, many of these models have been restricted to single 
isolated galaxies or cosmological zoom-in simulations of a 
very small sample of galaxies due to the expensive compu- 
tational cost of full cosmological simulations. 

Recently, both semi-analytic and non-equilibrium H2 
calculations have been implemented into full cosmological 
simulations. Kuhlcn et al. (2012) implemented the ana- 
lytic model of Krumholz et al. (2008, 2009) and McKee & 
Krumholz (2010) in the adaptive-mesh-refinement code Enzo 
(Bryan & Norman 1997; O'Shea et al. 2004) to study how 
Hb-based star formation affected dwarf galaxies at z > 4. 
Both their previous model and the new H2 model were able 
to reproduce many of the observational results pertaining to 
the KS relation. The advantage they found within the H2 
model was that it reduced the number of free parameters, 
and that star formation was quenched in dwarf galaxies from 
the onset without the need to artificially enhance stellar 
feedback. Christensen et al. (2012) implemented the non- 
equilibrium, fully time-dependent model of Gnedin et al. 
(2009) into their cosmological SPH code GASOLINE (Wad- 
sley et al. 2004) in order to study the effects of liVbased 
SF model on a dwarf galaxy down to z = 0. They found 
that the inclusion of H2 resulted in a greater baryonic mass 
in the disk, making it brighter, bluer, and more gas rich at 
z — than the same galaxy formed without the inclusion 
of H 2 . They also found that with H 2 there was more star 
formation at late times. 

While there are other models of star formation based on, 
for example, supersonic turbulence in the ISM (e.g. McKee 
& Ostriker 2007; Kritsuk & Norman 2011; Renaud et al. 
2012), it is still worthwhile to explore an implementation of 
H2-based SF as well, and investigate its implications. The 
purpose of this paper is not to decide which process triggers 
the star formation (i.e., supersonic turbulence or molecules), 
as our simulations have neither the resolution nor detailed 
dust physics to address the issue. In this paper, we limit 



ourselves to examining the effects of a new H2-based SF 
model on galaxy formation, and we defer the implementation 
of the turbulence-based SF model to the future. 

There is another good reason to study the H2-based 
SF model in cosmological simulations. Many of the earlier 
works based on a CDM model have predicted very steep 
faint-end of the mass/luminosity functions at high-redshift 
(e.g., Nagamine et al. 2004c; Night et al. 2006; Lo Faro 
et al. 2009; Finlator et al. 2011; Jaacks et al. 2012a), and 
suggested that these low-mass galaxies are responsible for 
reionizing the Universe at z ^ 6. However, the observational 
estimates yield slightly shallower faint-end slopes, and if the 
observational results are not affected by the magnitude limit 
very much the simulations need to consider processes that 
would decrease the number of low-mass galaxies, especially 
at high redshift. One of such candidate process is H2-based 
star formation, and Jaacks et al. (2012a) for example have 
speculated that the H2-based SF model may reduce the dis- 
crepancy in GSMF at the low-mass end. 

This paper is organized as follows. In Section 2 we dis- 
cuss simulation parameters, SF models, and basic results. 
Section 3 contains our findings for galaxy populations. The 
results of SHMR, cosmic SFRD, GSMF, and KS relation arc 
presented along with resolution studies. Lastly in Section 4 
we summarize our results and discuss future prospects. 



2 SIMULATIONS & BASIC RESULTS 

For our simulations we use a modified version of the 
GADGET-3 cosmological SPH code (originally described in 
Springcl 2005). Our conventional code includes radiative 
cooling by H, He, and metals (Choi & Nagamine 2009), 
heating by a uniform UV background (UVB) of a modified 
Haardt &; Madau (1996) spectrum (Katz et al. 1996; Dave 
et al. 1999; Faucher-Giguere et al. 2009), supernova (SN) 
feedback, the Multi-component Variable Velocity (MVV) 
wind model (Choi & Nagamine 2011), and a sub-resolution 
model of multiphase ISM (Springel & Hernquist 2003). In 
this multiphase ISM model, the high-density ISM is pictured 
to be a two-phase fluid consisting of cold clouds in pres- 
sure equilibrium with a hot ambient phase. Thermodynamic 
forces are calculated only for the hot phase. The cold phase 
on the other hand provides material for star formation, is 
subject to gravity, adds inertia, and participates in mass & 
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energy exchange with the hot phase at the sub-particle level. 
The primary purpose of this work is to improve upon the 
SF models implemented within this code. Previous SF model 
implementations are also discussed in upcoming sections. 

The same set of initial conditions (ICs) used by Choi & 
Nagamine (2011) and Jaacks et al. (2012a) are employed in 
this study to allow for cross comparison. Three primary sim- 
ulations were run consisting of 2x400 3 or 2x600 3 particles 
of gas and dark matter. Comoving box sizes of 10/i _1 Mpc, 
34/?T 1 Mpc, & 100ft _1 Mpc are used to capture a range of 
halo masses. These runs will be referred to by their par- 
ticle count followed by the length of their box: N400L10, 
N400L34, & N600L100. The other three runs (N144L10, 
N500L34, N600L10) where used mainly for resolution test 
of the H 2 model, and N500L34 & N600L10 runs were per- 
formed only for the H2 run due to expensive computational 
cost. The ICs were generated using cosmological parameters 
consistent with the WMAP best-fit values (Komatsu et al. 
2011): fi m = 0.26, Q A = 0.74, D. b = 0.044, H o /W0 = 0.72, 
n s — 0.96, and cr 8 = 0.80. The simulation parameter values 
are summarized in Table 1. The runs with smaller box sizes 
are stopped at earlier times, because they do not include 
longer wavelength perturbations. 

There are three additional differences between what we 
will refer to as the 'Fiducial' runs (Choi & Nagamine 2011; 
Jaacks et al. 2012a) and the 'H2' runs. First, we use an up- 
dated model of UVB radiation in the H2 runs, as will be dis- 
cussed in Section 2.2.4. The Fiducial run uses an older model 
wherein the UVB turns on at z — 6.08 to mimic the sudden 
reionization of the Universe; the UVB of the updated model 
turns on at z — 10.75, given the more recent WMAP results 
on early reionization of the Universe. For the majority of our 
comparisons, this change has little impact. Secondly, our H2 
runs use an optically-thick ultra-violet threshold or 'OTUV 
(Nagamine et al. 2010) which will be discussed in Section 2.3. 
Most comparisons presented in this paper are not effected 
by this; the column density distribution presented in Sec- 
tion 3.5 and Figure 18 however, is strongly effected by this 
change. Our low resolution (N144L10) Fiducial run uses the 
OTUV threshold. Lastly, the Fiducial run uses the Salpeter 
(1955) initial mass function (IMF), while our new runs use 
the Chabrier (2003) IMF. The choice of the IMF is reflected 
in the value of gas recycling fraction parameter /3 in the SF 
relation in our simulation. We have also verified that this 
has little impact on the results presented in this paper. 



2.1 Previous SF models 

2.1.1 SH model 

Springel & Hernquist (2003, SH model) describes the hy- 
brid multiphase model for SF, originally implemented in the 
GADGET code. In this model, the SFR is determined by 



(2) 



(1-/3) 



Pcold 



(1) 



where p co id is the density of cold gas available to form stars, 
and j3 is the fraction of stars instantaneously destroyed as 
supernova, determined by the stellar IMF. The parameter 
t, is the SF time-scale which is taken to be proportional to 
the local dynamical time of the gas: 



( P \ ~ 1/2 
i* = t I — , 

\PthJ 

where p th is a density threshold, above which the gas is as- 
sumed to develop a multiphase structure and form stars. 
The parameter to controls the normalization of the Kenni- 
cutt (1998) relation: 



(3) 



if S g 

jS fr-1 A ( Sgas/ iM pc- 2 )" if S g 

where T,th is the SF threshold surface density. Observations 
suggest A = 2.5 ± 0.7 x lO^Moyr^kpc- 2 , n = 1.4 ± 0.15, 
and E ih ~ lOM pc~ 2 (Kennicutt 1998). The cutoff in S SFR 
is controlled by p t h, which indirectly regulates T, t h- See 
Springel & Hernquist (2003) for a description of how p t h 
is calculated self-consistently for this model. 

Our simulations deal with three dimensional densities 
(i.e. Equation 1) rather than the two dimensional surface 
densities described by the KS-relation. It is not obvious then 
that Eq. (1) would be able to reproduce the observed KS- 
relation given by Eq. (3). Previous simulations (Springel & 
Hernquist 2003; Nagamine et al. 2004b) were able to demon- 
strate that the observed relation could indeed be reproduced 
using t ( * = 2.1 Gyr, which resulted in a threshold value of 
pth = 0.35 /i 2 cm~ 3 . However Nagamine et al. (2004b) and 
Choi & Nagamine (2010) showed that using this value of 
Pth results in overprediction of Ssfr at low column densi- 
ties (Nhi^ 10 20 ' 5 cm -2 ). This over-prediction, coupled with 
the fact that this model produces stars from cold atomic 
gas rather than molecular, leads to the necessity for im- 
provement in the sub-grid SF model, which we describe in 
the following sections. 

2.1.2 Pressure model 

Previous SF models assumed that the exponents of the Ken- 
nicutt and Schmidt relationships are equal. This is only true 
if the galactic disk scale- height is constant, or the equation 
of state behaves as P oc p 2 (Schaye & Dalla Vecchia 2008). 
Arguing that these assumptions are unnecessary and often 
incorrect, Schaye & Dalla Vecchia (2008) formulated a fully 
analytic conversion from the 2D KS-relation (E gas ) to a 3D 
Schmidt-relation (p ga s)- They proposed that the density of 
a self-gravitating disk will fluctuate on the local Jeans scale, 
leading to the scale-height also being on the order of the lo- 
cal Jeans scale. This in turn leads to the gas column density 
being on the order of the 'Jeans column density': 



^gas,,7 — Pgas-^J — 



7 V 7 %f p v/ 2 

q J (Jg^tot) , 



(4) 



where Lj = c s /^/Gptot is the Jeans length, c s — 
x/fPtot/ pgas is the local sound speed, 7 is the ratio of spe- 
cific heats, G is the gravitational constant, f g is the mass 
fraction within the scale-height of the gas, and Ptot is the 
mid-plane pressure. The SFR in this model is also described 
by Equation (1); the difference comes in the formulation of 
to which is derived by combining Equations (3) & (4) : 



tn = 



A' 1 (lM QP c- 2 ) n (lf s P tot ) 



(l-rO/2 



■. i • \ , , ■ ■ j (5) 

2jSFR V Or / 

Schaye & Dalla Vecchia (2008) claim that their analytical 
conversion renders their parameters 'untweakable'. Adopt- 
ing n = 1.4 & 7 = 5/3, Choi & Nagamine (2010) imple- 
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mented this model within our GADGET-3 code with some 
minor modifications. It was found to reduce the overpre- 
diction of Ssfr at low column densities and was in good 
agreement with the observed KS-relation. It also reduced 
the SFR in low-density regions, causing a suppression of 
early star formation, which in turn shifted the peak of the 
cosmic SFRD to lower redshifts in better agreement with 
observational estimates. 

The disadvantage of this model is that we are still im- 
posing the KS relation onto our SF prescription. In an ideal 
situation the model would naturally reproduce the observed 
KS relation without such impositions. Simulations and data 
from previous work based off of the Schaye & Dalla Vecchia 
(2008) Pressure SF model (Choi & Nagamine 2010, 2011; 
Jaacks et al. 2012a) will serve as our Fiducial runs for com- 
parison. 



2.2 H2 regulated star formation 

If star formation really requires molecular gas, then track- 
ing the H2 gas fraction and basing our SF prescription on 
it would make for a more realistic sub-grid model. There 
are currently two primary ways in which the H2 mass frac- 
tion can be tracked in simulations. The first is a non- 
equilibrium model which calculates the H2 fraction via a 
fully time-dependent chemistry and radiative transfer cal- 
culation as was done by Gnedin et al. (2009) and Chris- 
tensen et al. (2012). This approach accurately calculates the 
instantaneous H2 gas fraction, but can be computationally 
expensive. The second is a analytical approach developed by 
Krumholz et al. (2008, 2009) and McKee & Krumholz (2010) 
(hereafter KMT model), which calculates an equilibrium H2 
fraction assuming a formation-dissociation balance. 

These two methods were directly compared in 
Krumholz & Gnedin (2011); the second method was found to 
be a viable and nearly cost free alternative to the computa- 
tionally expensive first option at metallicities Z 10 _ Zq, 
where Zq is the solar metallicity. At metallicities < 10~ 2 Zq 
the KMT model was found to over-predict the fractional 
H2 abundance due to the neglect of time-dependent effects. 
Krumholz (2012) however, argues that the equilibrium H2 
fraction rather than the instantaneous one correlates more 
with gas temperature. He argued that the thermal timescale 
of gas is much shorter than the chemical timescale, which 
means that low metallicity clouds should cool via collisional 
de-excitation and form stars faster than they can fully con- 
vert all of their atomic gas to molecular. This suggests that 
the fractional H2 abundance calculated by the KMT model 
may more accurately predict the amount of material avail- 
able to form stars in low metallicity gas. Due to the compu- 
tational simplicity we choose to adopt the KMT model to 
track H2 within our simulations, and examine its impact on 
cosmological galaxy formation. 



2.2.1 KMT model 

Krumholz et al. (2008, 2009) and McKee & Krumholz (2010) 
developed an analytic model for tracking H2 mass frac- 
tion. They used a Stromgren-type analysis, starting with a 
spherical gas cloud immersed in a uniform, isotropic Lyman- 
Werner band radiation field. Assuming that the cloud is in 



a steady state, they proceeded to solve the radiative trans- 
fer and H2 formation-dissociation balance equations. After 
some approximations, they found a solution 



/h 2 



Ph 2 
Phi 



1 + 0.25s 



(6) 



where /h 2 is the H2 mass fraction relative to the neutral hy- 
drogen gas. At such high densities where molecular gas may 
form, the hydrogen gas would be mostly neutral within our 
multiphase ISM model, hence the neutral hydrogen (Hi) in 
the denominator above (see also Section 2.2.2). The param- 
eter s is essentially the size of the atomic-molecular complex 
given by 



tr 



ln(l + 0.6x + 0.01x 2 ) 
0.6r c 



(7) 



where tr is the dust optical-depth of the atomic-molecular 
complex, t c is the mean optical depth, and £d is the char- 
acteristic radius (in units of the cloud radius) at which the 
transition from dust-dominated to molecular-dominated ab- 
sorption occurs. The variable \ can be thought of as an 
estimation of the local radiation field given by 

°~d, -21 \ G' {) 
TZ-16.5 J fin 



x = n 



(8) 



Here cr di _2i is the dust cross section per H nucleus to 1000A 
radiation normalized to a value of 10 -21 cm 2 , 1Z-ie. 5 is the 
rate coefficient for H2 formation on dust grains normalized 
to the Milky Way value of 10 _16 5 cm" 3 s _1 (Wolfire et al. 
2008), G' is the ambient UV radiation field intensity nor- 
malized to the Draine (1978) value for the Milky Way, and 
tih is the volume density of H nuclei. 

At this point /h 2 depends only on r c & \- ^ n order to 
calculate the dust optical depth (r c ) we first calculate the 
local Hi column density (Shi) by means of a Sobolev-like 
approximation (e.g. Gnedin et al. 2009; Krumholz & Gnedin 
2011): 



Shi « Phi x h, 



(9) 



where h is the local scale height calculated by 

*"i^r (I0 » 

We find that per-particle, the Sobolev approximation gives 
higher values of Shi than a true column density calculation. 
However, when we take the mass- weighted average of these 
values along a column, the Sobolev approximation was ac- 
tually lower by a factor of ~ 5 in the high density regions 
of interest. This suggests that, within each column, there 
are more particles with a low Sobolev-surface density. The 
average therefore, is biased towards lower values. As an ex- 
periment we ran a test simulation boosting Shi in the high 
density regions by a factor of 5, and we find that more stars 
are formed at late times due to higher S values. Given that 
our current Sobolev approximation gives lower S values, we 
could regard our results on star formation as lower limits 
(See Section 3.2). 

We can then find the dust optical depth by the relation 
t c — Sgasfd/zxH, where ad is the dust cross section per hy- 
drogen nucleus and Ph is the mean mass per H nucleus. The 



dust cross section is taken to be o d _ 2 i 



10" 



where Z sn is the gas metallicity normalized to the solar 
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Figure 1. The parameter s (Equation 11) as a function of gas 
surface density for different metallicities. The value of s = 2 cor- 
responds to the transition from fully atomic gas to atomic & 
molecular within the KMT model. Lower metallicity gas requires 
larger column densities (i.e., more shielding) in order to form H2. 



neighborhood value Zq = 0.0204 (Rodriguez & Delgado- 
Inglada 2011). The KMT model shows that, if the ISM is in a 
self-consistently determined two-phase equilibrium, then the 
ratio G'o/nji takes a characteristic value, and Equations (7) 
& (8) become 



ln(l + 0.6x + 0.01 X 2 ) 



0.04 



M pc- 



and 



1 + 3.1(Z/Z Q )° 
4.1 



(11) 



(12) 



respectively. Using Equations (11) & (12) we can now calcu- 
late the H2 fraction of each gas particle by means of Equa- 
tion (6) if < s < 2 (McKee & Krumholz 2010), otherwise 
/h 2 = 0. 

Figure 1 shows how the transition from fully atomic to 
atomic & molecular changes for different metallicities. The 
value s = 2 is the minimum complex size required to al- 
low for the transition from atomic to molecular gas for any 
given metallicity. This represents the need for the gas cloud 
to be sufficiently large to allow for the formation of H2. 
External UV radiation is absorbed first by dust, which is 
essentially traced by the metallicity, then by a thin molec- 
ular region. If the cloud is too small then there will not be 
enough surrounding material to absorb all of the UV, and 
the Eb-core will be dissociated. If the cloud is large, we ex- 
pect a large molecular core. Higher metal content effectively 
absorbs more radiation, allowing for the formation of H2 at 
lower surface densities. 

Knowing /h 2 allows us to modify our SF model by 
replacing p co id in Equation (1) with the H2 mass density 
Ph 2 — /h 2 Phi- We also change our SF time-scale to a more 
physically realistic value, namely the free-fall-time of the H2 
gas available to form stars 



3tt 



32Gph 2 



(13) 



Furthermore, observations have shown that SF is a slow pro- 
cess and that the efficiency at which dense regions produce 
stars is ~1% (Krumholz & Tan 2007; Lada et al. 2010). 
To account for this we introduce an efficiency parameter of 
£SF = 0.01 which leads us to our new formulation of Equa- 
tion (1): 

Ph 2 



p, = (1 - 0) e SF 



If a gas particle has /h 2 
at the above rate. 



(14) 

> 0, then it is eligible to form stars 



2.2.2 Assumption on the neutral fraction 

As discussed in the previous section, in order to calculate the 
fractional H2 abundance, we must first calculate the scale- 
height of Hi, which then allows for the calculation of Em- 
Our GADGET-3 code tracks the neutral fraction of each gas 
particle individually. For the high density multiphase gas 
however, the neutral fraction is tracked only for the hot 
phase, and the cold gas fraction x = p c / p is computed within 
the multiphase ISM sub-particle model (Springel & Hern- 
quist 2003). For the very high-density particles, most of the 
mass is in cold, neutral phase (x > 0.95), but the tenu- 
ous hot phase determines the thermodynamic temperature. 
We calculate the neutral fraction using the x-parameter 
for high-density particles above the SF threshold for our 
N144L10 fiducial run at z = 3, and find all of particles to 
have /hi > 0.96, and 98% of the hydrogen mass to have 
/hi > 0.97. Given the small mass fraction of ionized gas, 
it is a good approximation to assume that any gas particle 
with nfjf > 0.6 cm -3 (Choi & Nagamine 2010) is completely 
neutral (/hi = 1) for the scale-height and column density 
calculations detailed in Section 2.2.1. 



2.2.3 H2 formation threshold density and Wind model 
modifications 

This new SF model (Eq. [14]) allows us to compute the den- 
sity threshold (pth) for individual particles based on their 
metallicity, above which results in the formation of H2. In 
the earlier version of our GADGET-3 code, Choi & Nagamine 
(2011) implemented the "Multi-component Variable Veloc- 
ity" wind model, in which a particle was allowed to travel 
as a wind particle with no hydrodynamic forces applied as 
long as its physical density exceeded n§f . The wind velocity 
of each particle was calculated based on the SFR of galaxy 
that the particle belongs to. We can now revise this wind 
criteria to be based off of individual particle's H2 formation 
threshold rather than a fixed value as in previous SF models. 

The formation of H2 requires sufficient shielding, or else 
the molecule will be dissociated by UV radiation. We can set 
the threshold for H2 formation for each particle by solving 
Equation (11) for E gas using s = 2; this value allows us to 
calculate the minimum E for SF within our model: 



ln(l + 0.6x + 0.01x 2 



(15) 



M pc- 2 2xO.O4(Z/Z ) 

We can then convert this surface density threshold to a 
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Figure 2. Physical density threshold (Eq. [16]) for H2 formation 
of all particles within a low resolution run (N144L10) at z = 3. 
Color gradient corresponds to the median scale-height h (Eq. [10]) 
as indicated by the color bar. The black dashed line represents 
the physical SF density threshold of nS? = 0.6 cm -3 used in our 
previous SF models (Choi & Nagamine 2010) . The KMT model 
generally requires higher densities to form H2 and hence be eligi- 
ble for SF, compared to the Fiducial model. 



volume density threshold using Equation (9) for each parti- 
cle, which leads us to the minimum volume density required 
to form H2 at that particle's particular metallicity: 



In 



Phi,- 



(l + 0.6x + 0.01x 2 )M Q pc- 
2 x 0.04 (Z/Z e )h 



(16) 



In other words, this is the H2 formation density threshold. In 
the present work, if the physical density of a gas particle is 
above its particular adaptive H2 formation threshold pm th , 
then it is eligible to be a member of the wind. 

Figure 2 shows the values of PHi th vs. metallicity of 
each particle in a low resolution simulation (N144L10). The 
previous SF density threshold is shown as the black dashed 
vertical line. Generally, the values of pm th are higher for 
each particle than in the previous SF models, allowing for 
particles to reach higher densities before becoming eligible 
to form stars. This makes SF in the H2 model less efficient 
than in the previous SF models. This plot also shows that for 
a given metallicity, a lower h returns higher pHi th • Particles 
with higher metal content have lower formation thresholds 
due to their ability to absorb more dissociating photons. 
The accumulation of particles at Z = 1O _3 Z0 corresponds 
to those that have yet to be enriched by SF, but they have 
varying pHi th due to variations in h. Some particles at — 2 < 
log n t h < have already been enriched by z — 3. 



2.2.4 Metal Floor 

In our Fiducial runs, initially all gas particles are metal free. 
Enrichment occurs during star formation; in this process SN 
explosions return a metal mass of AMz = J/.AM, to the 
ISM, where y„ — 0.02 is the yield, and M* is the mass of 
newly formed long-lived stars. It is assumed that each gas 
particle behaves as a closed box locally, wherein metals are 



instantaneously mixed between cold clouds and ambient hot 
gas. 

Within the framework of our new SF model, stars can 
only form if they contain H 2 , as determined by Equation (6). 
In order for fn 2 7^ 0, the gas particle must have some metal 
content. To begin forming stars, we must first enrich the gas 
particles by hand at very high redshift. Recent high reso- 
lution numerical studies by Wise et al. (2012) found that a 
single pair-instability supernova of a Pop. Ill star can enrich 
its host halo to a metallicity of 1O _3 Z0. Their findings are 
in agreement with observed DLA metallicities, which have 
metal floors of the same order (Wolfe et al. 2005; Penprase 
et al. 2010). To mimic this enrichment, we introduce a metal 
floor of Z = 1O _3 Z0 for all gas particles at a specified epoch. 

To test the impact of the assumed metal floor, a few low 
resolution simulations (N144L10) are performed introduc- 
ing the metal floor at redshifts of z = 9, 11, &13; we refer to 
these as MF9, MF11, & MF13 runs, respectively. The cosmic 
SFRD histories are nearly identical between these three sim- 
ulations; they differ only in the point at which each simula- 
tion starts to form stars. This is due to the different times at 
which their metal floors are introduced. The MF11 & MF13 
runs both start to form stars at z ~ 9.2, while MF9 does not 
begin star formation until z ~ 8.6. After their initial bursts 
of star formation, each of the three simulations begin follow- 
ing the same SFRD track from z ~ 8.3 to z = 3. The GSMF 
between the three simulations are also nearly identical at 
z — 3 & 6, suggesting that the formation of galaxies within 
our simulations does not heavily depend on when the metal 
floor is set. Lastly, the SHMR (cf. Section 3.1.2) is examined 
at z = 3 & 6 for each simulation. We find significant scat- 
ter in the results for all three runs, but the median SHMR 
values for each simulation are all well within one standard 
deviation of one another. This again suggests that the stellar 
fraction of each halo does not depend heavily on the time at 
which the metal floor is set. 

The redshift at which the metal floor is introduced is 
related to the metal enrichment by Pop. Ill stars. There- 
fore, we choose to introduce our metal floor at the epoch 
of reionization. Observational estimates by Komatsu et al. 
(2011) suggest this happens at redshift z — 10.6 ± 1.2. In 
our simulation, reionization is set by the UV background 
model (Faucher-Giguere et al. 2009, December 2011 ver- 
sion) 1 , hence our metal floor of Z = 10~ 3 Z Q is set at 
z ~ 10.75 accordingly. 



2.3 Gas phase diagram 

We examine a low resolution N144L10 simulation to study 
the gas temperature-density phase diagram. Figure 3 com- 
pares our Fiducial run with the new H2 run at z = 3. 
Panel (a) represents our Fiducial run and contains three ver- 
tical lines representing different thresholds. The left most 
dashed line (n^ h v = 0.006 cm 3 ) represents an optically- 
thick ultra-violet threshold or 'OTUV (Nagamine et al. 
2010). Any particle below this threshold will be heated by 
the UVB to > 10 4 K; any particle above this threshold is 
shielded from the UVB. Nagamine et al. (2010) and Yajima 



https : //www . cf a. harvard . edu/~cgiguere/UVB .html 
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Figure 3. Gas temperature vs. number density phase diagrams for the low resolution N144L10 runs at z = 3. Panel (a): the Fiducial 
run. The right-most dotted line is the physical SF and wind density threshold (n§F = 0.6cm -3 , Section 2.3), the middle dot-dashed line 
is the maximum wind travel length (TL) discussed in Section 2.3, and the left-most dashed line corresponds to the OTUV threshold 
discussed in Section 2.3. Panel (b): H2 run. Here only the OTUV threshold is shown. There is no fixed SF density threshold, as it is 
different for every particle depending on the metallicity and surface density as described in Section 2.2.3. 



et al. (2011) have demonstrated that this threshold is rea- 
sonable using full radiative transfer calculations. 

The right most dotted line in Figure 3a represents the 
constant SF physical density threshold (nfjf = 0.6 cm -3 ) in 
the Fiducial run. Any particle whose density exceeds this 
threshold is allowed to form stars based on the prescrip- 
tion outlined in Section 2.1.2. At densities & temperatures 
above n ~ 3 cm' " 3 & T ~ 10 4 K, we begin to see the effects 
of SN feedback and the multi-phase ISM model (Springel & 
Hernquist 2003). The cold phase component dominates the 
mass of the particle, but the hot phase governs the temper- 
ature. What is seen in this 'arm' is the temperature of the 
gas heated by SN feedback (hot phase component), and the 
density of the cold phase component. When we direct our 
attention to the H2 run (Panel b) we see that this arm is 
now extended out to higher densities at lower temperatures. 
This is because stars are only allowed to form if the gas par- 
ticles contain any H2 above the adaptive density threshold 
Pm,th given by Equation (16). As discussed earlier, Figure 2 
illustrates how pm.th is typically higher than the previous 
SF density threshold, which allows particles to condense to 
higher densities without being heated by SN feedback. 

The dot-dashed line in between the two previously dis- 
cussed lines in Figure 3a represents the maximum wind 
travel length (TL) threshold of n^ L = 0.1xn t s ,f = 0.06 cm -3 . 
If a particle becomes a member of the wind, hydrodynamic 
forces are turned off for a brief period of time (Springel & 
Hernquist 2003; Choi & Nagamine 2011). This allows the 
gas to adiabatically expand and cool to lower temperatures 
until the density drops below h"l, or the brief period of time 
has elapsed. The dot-dashed line is absent from the Panel 
(b) because of the varying density threshold for each parti- 
cle. Instead of a temperature discontinuity, as can be seen 




20 21 22 20 21 22 20 21 22 
Log S ffl +H, [cm -2 ] 

Figure 4. H2 mass fraction as a function of H1+H2 surface den- 
sity within our N600L10 run (black points) at z = for three dif- 
ferent metallicity ranges. The red shaded regions show the results 
of three Milky Way-like simulations of Christensen et al. (2012) 
using a full non-equilibrium H2 model with different metallicities 
of (a) 1Z Q , (b) O.3Z , and (c) 0.1Z e . 



in Panel (a), we have a 'tail' that extends all the way to 
the OTUV threshold. This tail corresponds to wind parti- 
cles that have adiabatically expanded as part of the galactic 
wind. 
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2.4 Atomic to molecular transition 

It is important to study where the atomic to molecular tran- 
sition occurs in our simulations. Figure 4 shows this tran- 
sition as a function of gas surface density in our N600L100 
run at z — for three different metallicity ranges. In the 
KMT model, the value of /h 2 is solely determined by the 
surface density of gas and metallicity (Equations 6, 11 & 
12), therefore the simulation data (black dots) in all three 
panels is restricted to a relatively thin band determined by 
the range of metallicity chosen. 

Christensen et al. (2012) examined this transition for 
three isolated Milky Way-like simulations at different metal- 
licities to test their fully time-dependent, non-equilibrium 
H2 calculation. Their raw simulation output can be seen as 
the red contour in Figure 4. The transition in our simula- 
tions is in good agreement with theirs, corroborating the 
comparison work of Krumholz & Gnedin (2011). However, 
our data deviates to higher column densities at high molec- 
ular fractions due to the per-particle overestimation of the 
column density calculated by the Sobolev-like approxima- 
tion, as discussed in Section 2.2. This deviation should not 
pose a problem since particles in these high column density 
regions are already fully molecular. 



3 RESULTS ON GALAXY POPULATIONS 
3.1 Dark matter halo content 

Dark matter (DM) particles were grouped using a simplified 
version of the parallel friends-of-friends group finder SUB- 
FIND (Springel et al. 2001). The code groups the particles 
into DM halos if they lie within a specified linking length. We 
adopt a standard value of b — 0.2 for the linking length pa- 
rameter, which is a fraction of the initial mean inter-particle 
separation. Each halo must also have a minimum of 100 par- 
ticles to be considered a halo. 

The DM halo mass functions were found to be in agree- 
ment between the H2 and Fiducial runs. This is an ex- 
pected result, because both sets of simulations were started 
from the identical IC. Both results are in good agreement 
with the Sheth & Tormen (1999) mass function. This paper 
focuses on baryonic properties, and it is useful to exam- 
ine and compare the contents of these halos between the 
Fiducial and H2 runs. The contents of each halo are calcu- 
lated by the summation of particle properties located within 
r 2 oo = [(GM DM )/(!!4z)J(z) 2 )] 1/3 (Mo & White 2002) of 
each halo's center of mass. 

3.1.1 Baryon fraction 

Figure 5 presents the baryon mass fraction over halo mass 
(/b = M b aryon/M D M = (Mgas,200 + M*, 2 oo)/Mdm) as a func- 
tion of log Mdm at z = 6,3,1 & 0. Here the cosmic mean 
(fib /Odm) is illustrated by the horizontal black dashed line. 
Panel (a) shows the composite data from the N400L10, 
N400L34, & N600L100 runs at z = 6; Panel (b) is com- 
posite data from the N400L34 & N600L100 runs at z = 3; 
Panels (c) & (d) are composed of data from the N600L100 
simulation only at z — 1 and 0, respectively. Solid lines rep- 
resent the median value within each bin, and the hatched 
regions represent a la spread in the data. The cutoff of the 
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Figure 5. Baryon mass fraction within r2oo of each halo, 
fb = M ha , Tyon /M DM = (M gas ,200 + M*,20o)/Mdm, as a function 
of halo mass (only DM) for 2 = 6, 3, 1, & 0. The red and blue 
solid lines represent the median points in each mass bin for the 
Fiducial and H2 runs, respectively. The hatched regions represent 
1(7 scatter in each Mdm bin. The cosmic mean baryonic fraction 
(^b/^DM) i s represented by the dashed horizontal line. 

data at lower mass end is determined by the halo mass limit 
of the halo grouping. 

At z — 6 (panel [a]), the /b of the two SF models agree 
with each other well and with the cosmic mean for halos 
above Mdm ~ 10 9 Mq. Halos below this mass in the H2 run 
have lower /b than in the Fiducial run by ~ 35%. This is 
presumably due to the different UVB model between the 
two runs; in the H2 run the UVB is turned on much earlier, 
and the gas in low-mass halos are photo-evaporated. 

This large difference in /b is nonexistent in low mass 
halos at z — 3 as shown in Panel (b). The median values 
within the H2 run are now generally higher than those in 
the Fiducial run. As we will see in later sections, this is 
likely due to higher SFRs and hence stronger SN feedback 
in the Fiducial run, and this trend continues to z = 0. The 
scatter in /b at Mdm ~ 10 9 ' 7 Mq is greater for the H2 model, 
but it encompasses that of Fiducial run. Both begin to drop 
slightly below the cosmic mean at around Mdm ~ 10 12 Mq. 

By z = 1, /b in the most massive halos settle to a 
value that is lower than the cosmic mean by ~ 7%. Again in 
Panel (c), we see /b in massive halos with Mdm > 10 12 M Q 
is in agreement between the two models. At lower Mdm, the 
Fiducial run still shows a lower baryon fraction. Finally at 
z — 0, we see that, for the halos with Mdm < 10 12 Mq, the 
mean /b has decreased slightly since z = 1. This means that 
the halos on average have acquired more dark matter than 
baryons (either through mergers or accretion), and/or at the 
same time lose the gas from galaxies due to SN feedback. The 
scatter of fb for lower mass halos is generally greater than 
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Figure 6. The SHMR as a function of total halo mass (DM+baryons) within T2oo- The data from semi-analytic models and observations 
are shown as the grey shade, which is identical in all four panels as it doesn't evolve very much with redshift (Behroozi et al. 2012). 



for the massive halos, and this is a known resolution effect 
(e.g., O'Shea et al. 2005). 



3.1.2 Stellar-to-halo mass ratio (SHMR) 

The ratio of stellar-to-halo mass as a function of total halo 
mass Mtot,200 (DM+gas+stars within V2oo) provides a use- 
ful insight on the efficiency of star formation in different 
halos, and it has collected significant attention in the re- 
cent years (Conroy et al. 2007; Baldry et al. 2008; Behroozi 
et al. 2010; Moster et al. 2010; Foucaud et al. 2010; Evoli 
et al. 2011; Leauthaud et al. 2012; Papastergis et al. 2012). 
All of these work find a prominent peak in this relation at 
Mtot,200 ~ 10 12 Mq, suggesting that there is a characteristic 
halo mass that galaxy formation is most efficient. This mass- 
scale roughly corresponds to the characteristic stellar mass 
M* of GSMF, i.e., the knee of Schechter function, where 
most of the stellar mass has formed globally. To further sur- 
prise, Behroozi et al. (2012, hereafter B12) found that, us- 
ing observational data and semi-analytic modeling, SHMR 
is almost time-independent between z = 4 to z = 0. This 
is interesting, because SHMR should reflect all the cumu- 
lative effects of past star formation and feedback, and non- 
evolving SHMR suggests tight self-regulation and a subtle 



balance between star formation and feedback. It is a signif- 
icant challenge for any cosmological hydrodynamic simula- 
tion to reproduce this relation across a wide range of halo 
mass and cosmic time. 

Note that we are specifically using M to t,200, and not 
Mdm for this comparison. This is because Munshi et al. 
(2012) pointed out that, in order to accurately compare sim- 
ulations to semi-analytic model results (such as abundance 
matching technique), one must observe the simulations in 
a similar manner. They refer to the work by Sawala et al. 
(2012), who showed that M to t,200 in iV-body simulations can 
be greater than those in iV-body+hydro simulations by up 
to 30%, because various baryonic processes (gas pressure, 
reionisation, supernova feedback, stripping, and truncated 
accretion) can remove baryons from the halo, decreasing 
the total mass systematically. Additionally, Graham et al. 
(2005) found that the stellar mass component of observed 
galaxies could be systematically underestimated by ~ 20% 
in some cases; for example, additional flux of 0.22 mag lies 
beyond the SDSS Petrosian aperture for a galaxy that has 
an i? 1/4 profile. Here, we consider that it would be more 
natural to examine SHMR as a function of M to t,200 rather 
than correcting our results by a certain number. 

In Figure 6, we show the SHMR in our simulations by 



© 0000 RAS, MNRAS 000, 000-000 



10 Thompson, Nagamine, Jaacks, & Choi 




Figure 7. Cosmic SFRD for our simulations compared to some observations. The left panel is for the Salpeter IMF, and the right for 
Chabrier IMF. The Fiducial runs are using Salpeter IMF, and the H2 runs are using Chabrier IMF. The observational data are from: the 
CLASH program (Postman et al. 2012; Coe et al. 2012, red triangles), Bouwens et al. (2011, 2012, black circles), Reddy & Steidel (2009, 
green crosses), Schiminovich et al. (2005, black stars), Kistler et al. (2009, cyan shade), and Nagamine et al. (2006, magenta hatched 
region). All observational data are corrected for dust extinction by each author as they deemed appropriate. 



calculating the total stellar mass contained within r2oo of 
each halo's center-of-mass (M*,20o)- If we assume that the 
B12 data extends out to z = 6, we see in Figure 6a that 
our H2 run does a good job at reproducing the B12 data at 
Mtot,200 < 10 12 M©, much better than the Fiducial model. 
We have verified that the different UVB models do not im- 
pact this result. The comparison of the models in Figure 6 
clearly suggests that the suppression of star formation in 
low-mass halos is favorably achieved by the H2 model. Note 
that this SF suppression is not due to the SN feedback, 
but rather due to the metallicity dependence of the new 
Eb-based SF model. This could become a critical point to 
distinguish between the Eb-based and turbulence-based SF 
models in the future. 



At 2 = 3 (Figure 6b) the SHMR increases slightly for 
both simulations. Our simulations are in agreement with 
the B12 data at Mtot.,200 < 10 12,2 Mq. However, we con- 
tinue to overestimate SHMR at M to t,2oo > 1O 12 ' 2 M0 down 
to z = 0, which could be due to lack of AGN feedback in 
our current simulations. It is widely believed that both ther- 
mal and momentum feedback from supermassive black holes 
suppresses the star formation in massive halos, making them 
'red & dead' (e.g., Di Matteo et al. 2005; Springel et al. 2005; 
Nagamine et al. 2005; Ostriker et al. 2010; Choi et al. 2012). 
There is little evolution between z = 1 & in our simula- 
tions (Figure 6c, d), and our results are higher than B12 data 
even for low-mass halos at z ^ 1. 



3.2 Quantities related to star formation 

3.2.1 Cosmic star formation history 

With the H2 model producing less stars in lower mass ha- 
los, we expect the cosmic SFRD to be lower as well. When 
comparing simulations to observational estimates of SFRD, 
we have to be careful about which IMF is being used. The 
Fiducial and H2 runs use different IMFs. In order to fairly 
compare the two, we must make corrections to either the 
simulation data or the observations. Because SFR is a raw 
output of our simulation, we prefer to take the latter route 
and correct the observational data to the same IMF as in 
simulation. A simple factor /imp allows for this conversion: 

• IMF -Salpeter i p /i*r\ 
P* =P* //IMF, (17) 

where pi MF represent an arbitrary IMF. To convert from 
Salpeter to Chabrier, we take /imf = 1-6 (e.g., Nagamine 
et al. 2006; Raue & Meyer 2012), and from Salpeter to 
Kroupa we take /imf = 1-8 (Horiuchi et al. 2009). This 
is because, for a given amount of observed rest-frame UV 
flux, the required SFR would be lower for an IMF that is 
weighted more towards higher masses. After correcting our 
IMFs, we find that both simulations roughly agree with the 
observed data. 

Figure 7 shows the cosmic SFRD history as a func- 
tion of redshift. As expected, the H2 runs show significantly 
lower SFRD at most redshifts relative to the corresponding 
Fiducial runs. The H2 runs do not start forming stars un- 
til z ~ 10.5, which is a consequence of our model. In the 
H2 run, in order for gas to be eligible for SF, it must first 
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Figure 8. Relationship between the masses of simulated galaxies 
(identified by the friends-of-friends grouping) and their nearest 
DM halos. Note that M+ is not exactly same as We see 

that the low-mass galaxies with M* ~ 10 6 — lO 9 M0at 2 = 6 
reside in more massive DM halos in the H2 runs than in the 
Fiducial run. The dashed line in each panel represent the scaling 
of logAfrjM = 0.8(logM* — 10) + 12, which is a simple eye-ball 
fit at z = 0. 



contain H2, which requires non-zero metal content. As dis- 
cussed in Section 2.2.4, we introduce the metal floor by hand 
at z ~ 10.75, after which stars are able to form. Until the 
metal floor is introduced, the gas continues to condense to 
higher densities. 

The slopes of the H2 SFRDs at high-z are slightly 
steeper than the Fiducial runs, because the H2 run starts 
SF a little later than the Fiducial run, and has more abun- 
dant high-density gas available for SF. It tries to catch up 
to the Fiducial run after the metal floor is introduced. For 
the same reason, the peak of the SFRD of the N600L100 H2 
run is slightly shifted towards lower redshift compared to 
the Fiducial run. With a lower SFRD in the H2 runs, there 
is more gas available for SF at later times. 

Note that Figure 7 is only showing the results of differ- 
ent simulations as separate lines, and it does not show the 
expected total SFRD. To really obtain the total SFRD, we 
must carefully examine the contribution from each galaxy 
mass ranges to SFRD, and sum up the contribution from 
each simulation. This was done by Choi & Nagamine (2010) 
for the Fiducial runs, and we will present such analyses sep- 
arately (Jaacks et al. in preparation). Here, we simply note 
that the expected total SFRD would be even higher than the 
red dot-dashed line of N400L10 run, and it would go right 
through the data range shown by the cyan and magenta 
shaded regions. 



Figure 9. Gas mass fraction / gas = Af gas / (Af gas + M*) of simu- 
lated galaxies as a function of galaxy stellar mass. Black triangles 
in Panels (b) & (c) are observed galaxies at 2 ~ 2 (Erb ct al. 
2006). In Panel (d) we show observational data at 2 = taken 
from Peeples k, Shankar (2011). 

3.2.2 In which halos do galaxies sit? 

So far, we have not considered the grouping of galaxies them- 
selves (i.e., star and gas particles). To examine individual 
galaxies in our simulations, we group gas and star particles 
based on the baryonic density field rather than the dark 
matter. This allows us to identify galaxies directly, and then 
calculate properties such as their SFRs, stellar masses {M t , 
which we distinguish from M^oo), gas masses (M gas , which 
we distinguish from Af gas ,2oo) and metallicities. A more de- 
tailed description of this galaxy group finding process can 
be found in Nagamine et al. (2004c). 

Together with the friends-of-friends halo finding, we are 
interested in how the grouped galaxies relate to the DM 
halos. To find out the matching between the two sets, we 
search for the nearest DM halo from the center-of-mass of 
each galaxy. We limit our galaxy sample to those with mini- 
mum 10 star particles, and those that reside in halos with at 
least 100 DM particles. Note that the DM structure between 
the Fiducial and H2 runs are nearly identical, because they 
both use identical ICs. 

We can then make a scatter plot of corresponding Mdm 
and M* of each halo as shown in Figure 8. In Panel (a) we see 
that the low mass galaxies (M* ~ 10 6 — 10 9 M Q ) at z = 6 in 
the two runs reside in different halos with different masses; 
the median result of the two runs lie almost an order of 
magnitude apart, with the Fiducial run galaxies residing in 
lower mass halos on average. This is because the SF requires 
a higher threshold density in the H2 run and the gas requires 
a deeper potential-well of massive halos in order to form the 
same amount of stars as in the Fiducial run. The results of 
the two runs converge at higher masses (Mdm > 10 9 M Q ), 
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Figure 10. Gas metallicity of simulated galaxies as a function 
of galaxy stellar mass. Black dashed line in all panels is the the- 
oretical model from Savaglio et al. (2005). Colored contours in 
Panel (d) correspond to observational data from Tremonti et al. 
(2004, dark-gray), Kewley & Ellison (2008, magenta), and Lara- 
Lopez et al. (2012, yellow). 

suggesting that those halos contain similar galaxies in the 
two runs. For higher mass halos with Mdm > 10 11,5 Mq, 
there seems to be little difference in the SF between the two 
runs. 

This M* — Mdm relation does not evolve significantly as 
time proceeds. For comparison, the dashed lines in Figure 8 
represent the same scaling of log Mdm = 0.8(log M* — 10) + 
12. The figure shows that the halos grow in mass with time, 
and the median lines slide up to top-right direction along 
the dashed line. 



3.2.3 Gas mass fraction of simulated galaxies 

Figure 9 shows the median gas fraction (/ gas = 
M gas /(M gas + M*)) of simulated galaxies as a function of 
galaxy stellar mass. In general, / gas in the Lh run is higher 
than that in the Fiducial run, but the la regions overlap 
with each other. The non-smoothness in the median lines in 
Panels (a) & (b) are simply due to the procedure of combin- 
ing the data from multiple simulations with different reso- 
lution. 

We find that the median lines do not evolve very much 
over time. At 2 ^ 3, / gaa declines steeply with M*; from 
values close to unity at M* ~ 1O 9 M to / gas < 0.1 at M* ~ 
1O 12 M0. This suggests that the massive galaxies with M* ~ 
10 12 M Q in our simulations have converted most of baryons 
into stars, and not much gas is left in them, coinciding with 
the downturn of the SFRD at these epochs (Figure 7). 

Black triangles in Panels (b) & (c) are from a sample of 
galaxies at 2 ~ 2 (Erb et al. 2006). Simulated galaxies from 



the H2 run tend to agree better with the observed data at 
2 = 3 & 1. In Panel (d) we show observational data of nearby 
galaxies from McGaugh (2005, stars), Leroy et al. (2008, 
filled circles), and West et al. (2009, 2010, crosses). Neither 
the H2 or Fiducial models agree well with observations at 
2 = 0. This may in part be due to the limited mass resolution 
of the N600L100 run; a higher resolution run would resolve 
lower mass galaxies, possibly shifting the distribution to the 
left in better agreement with observations, if we were to 
make a composite plot from different runs. Another possible 
cause for this discrepancy is that too much unenriched gas 
has fallen into these massive galaxies between z — 1 and 
2 = 0, pushing /gas to higher values. 

3.2.4 Metallicity of galaxies 

As the gas recycling with metals takes place following star 
formation, the metallicity of galaxies should be roughly in- 
versely proportional to their gas fraction if we neglect gas 
infall. In Figure 10, we compute the average metallicity of 
each galaxy by summing the SFR-weighted metallicity of 
all gas particles within grouped galaxies. Observational con- 
straints on galaxy metallicities are derived from Hll regions, 
therefore only the regions close to bright stars are probed, 
and a SFR-weighted metallicity is more appropriate for com- 
parison than a mass-weighted metallicity. If we instead use 
a mass-weighted metallicity, we obtain much lower metal- 
licity than the SFR-weighted one, because it would include 
unenriched gas in the outskirts of galaxies. 

Figure 10 illustrates a general agreement with the above 
expectation: higher mass galaxies have lower / gas and higher 
metallicity in general. In all four panels, we show model pre- 
dictions from Savaglio et al. (2005) as black dotted lines. 
At 2 = 6 we find our simulations to be in rough agree- 
ment with the model at M* > 10 8 Mq, but below this mass 
we over-predict Savaglio's model result. At lower redshifts, 
the metallicity of simulated galaxies is always below that of 
Savaglio's model. The colored regions at 2 = correspond 
to observations from Tremonti et al. (2004, dark grey), Kew- 
ley & Ellison (2008, magenta), and Lara-Lopez et al. (2012, 
yellow). In particular, Kewley & Ellison (2008) have shown 
that the observational estimates of metallicity could vary 
significantly depending on the adopted estimator. The me- 
dian of our Fiducial run overlaps with the magenta shade in 
a limited stellar mass range while the scatter overlaps the 
entire range. Our H2 simulations however, are in excellent 
agreement with the magenta shade while under-predicting 
the metallicity of high mass galaxies with M* xi 1O 1O M0 
when compared to the dark grey and yellow shades. 

If we do not weight by SFR, we find that the mass- 
weighted metallicity is much lower for the H2 run, but the 
Fiducial run results are not so much affected. This is be- 
cause, for a fixed baryonic density cut used in our galaxy 
grouping procedure (Nagamine et al. 2004c), the galaxies in 
the H2 run will contain less star-forming gas than in the 
Fiducial run due to the higher Eb-formation (and hence SF) 
threshold, as we discussed in Figure 2. 

3.2.5 Specific star formation rates of galaxies 

Figure 11 shows the redshift evolution of the specific 
star formation rates (i.e., SFR per unit stellar mass, 
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Figure 11. Specific star formation rate (sSFR = SFR/M*) of simulated galaxies as a function of galaxy stellar mass. The observed 
data ranges are indicated by the shaded region with the source indicated in each panel. Contoured observational data was taken from 
Krumholz & Dekel (2011), while the observations at z = 6 were taken from Ono et al. (2010, circles) and Labbe et al. (2010b, triangles). 



sSFR = SFR/M*) in our simulations. This plot shows the 
instantaneous efficiency of SF, whereas the SHMR reflects 
all past history of SF and feedback. Panel (a) shows that the 
low mass galaxies in the H2 run at z — 6 have higher sSFRs 
than those in the Fiducial runs, mirroring the steeper slope 
of the SFRD (Figure 7) for the H2 run. Our simulation data 
is higher than the observational data of Lyman-a emitters 
at z = 5.7 & 6.6 (Ono et al. 2010) and z-dropout galaxies 
at z ~ 7 (Labbe et al. 2010b), but within their error bars. 

The H2 run in Panel (b) (z = 3) again show a slightly 
higher sSFR than the Fiducial run for lower mass galaxies 
with M* < 1O 9 ' 6 M0. At higher masses, the two runs agree 
very well, as well as with the observed data at z = 3.7 & 
2.0, indicated by the shaded region (Daddi et al. 2007; Lee 
et al. 2011). Panel (c) (z — 1) also shows similarly good 
agreement between the two runs and the observational data 
range (Elbaz et al. 2007). 

Panel (d) (z = 0) shows that the sSFR of both runs 
continue to decrease with time, but the H2 run decreases at 
a faster rate. Therefore the Fiducial run has a higher sSFR 
than the H2 model at z = 0. Both models agree with the 
observational data (Brinchmann et al. 2004) with a slightly 
decreasing sSFRD as a function of M*. 



3.2.6 Redshift evolution of the sSFR 

Observations indicate that galaxies of similar mass (~ 
1O 1O M0) have relatively constant sSFRs on the order of 
1 — 2Gyr _1 in the redshift range of 2 = 2 — 7 (e.g., Stark 
et al. 2009; Gonzalez et al. 2010; Labbe et al. 2010b,a). This 
sSFR 'plateau' is difficult to produce with current models of 
galaxy formation (e.g., Bouche et al. 2010; Weinmann et al. 
2011). Figure 12 shows the sSFR as a function of redshift 
for simulated galaxies with M* = 1O 1O M0. We find that the 
sSFR for these galaxies decline gradually rather than a steep 
drop off around z ~ 1 as observations suggest. The H2 run 
however, produces a slightly steeper drop off at z < 1 than 
the Fiducial run, but the normalization is still lower than 
the compilation of observed data points by Weinmann et al. 
(2011). Neither model produces the observed plateau in the 
redshift range of z — 2 — 7. 

For comparison, simulation data from two different 
wind models of Dave et al. (2011, dashed & dot dashed lines) 
are included, and they show similar trends to our simulation 
data. Our results are also very similar to the results of the 
semi-analytic model of Neistein & Weinmann (2010); Wein- 
mann et al. (2011) without the plateau at z > 2. The general 



© 0000 RAS, MNRAS 000, 000-000 



14 



Thompson, Nagamine, Jaacks, & Choi 



1.5 

1.0 

" 0.5 

! 

u_ 

10 -0.5 
i/i 

Dl 
O 



■1.0 
-1.5 
■2.0 










I f 















H2 

Fiducial 

Reddy et al. 2009 
Stark et al. 2012 
Weinmann et al. 2011 
Dave et al. 2011 SW 
Dave et al. 2011 VZW 



1 



Figure 12. Redshift evolution of the sSFR of simulated galaxies. 
Data points for the Fiducial and H2 runs are the median sSFR 
at Af* = 10 10 M Q (Figure 11), while the error bars represent a la- 
spread in the data. Observations are taken from Reddy & Stei- 
dcl (2009, stars), Weinmann et al. (2011, cyan shade), Bouwens 
et al. (2012, squares), and Stark et al. (2012, circles). Simulation 
data from Dave et al. (2011) is shown as the black dashed (VZW 
model) and dot-dashed (SW model) lines for comparison. Fiducial 
points are offset by 0.1 dex for clarity. 



agreement between multiple different simulations and semi- 
analytic models of galaxy formation suggest that the ACDM 
model predict a general decline in the sSFR of galaxies of 
a given mass, contrary to observations. However, we note 
that none of these simulations included the effect of AGN 
feedback. Krumholz & Dekel (2011) argued that taking the 
metallicity-dependence of H2 formation would help to recon- 
cile the discrepancy, however, even with our new H2-based 
SF model, our simulations do not produce the plateau of 
sSFR at 2 > 2. 



3.3 Galaxy stellar mass function (GSMF) 

In the previous sections, we have seen that SF is less efficient 
in the H2 run, which should also be reflected in the GSMF. 
Recall that for a given M* at high-z, the galaxies reside 
in more massive halos in the H2 runs (Figure 8). Since the 
higher mass halos are less abundant in a CDM universe, this 
will reduce the number of low-mass galaxies and shifts the 
galaxy population to higher mass DM halos. 

Figure 13 shows the GSMF for our three primary runs 
(N400L10, N400L34, N600L100) at 2 = 6. In Panels (a-c) we 
directly compare the H2 run to the corresponding Fiducial 
run for each simulation, and find that the H2 run produces 
far fewer low-mass galaxies as expected. Note the different 
y-axis ranges in Panels (a-c). Our result is in general agree- 
ment with the findings of Kuhlen et al. (2012); they also 
found a decrease in their GSMF at M* < 1O 9 M at z = 4. 

Figure 13d shows the comparison of the composite 
GSMF from the two runs, following the method of Jaacks 
et al. (2012a); we connect the GSMF from runs with differ- 
ent box sizes at the resolution limit of each run. This method 
allows us to cover a wider range of M+ utilizing many simula- 
tions, and present the results collectively. The observational 



estimate from Gonzalez et al. (2010, yellow shade) at 2 = 6 
is also shown. At the high-mass end of M* > 10 9 Mq, the 
two composite GSMFs from H2 and Fiducial runs agree well. 
The slight kink in the composite GSMF at M* ~ 10 8 ' 8 M Q 
for the H2 run is due to the resolution gap between the sim- 
ulations; we have verified that an intermediate resolution 
run (N500L34, e = 2.72 ft _1 kpc) fills in this gap. Due to the 
heavy computational load, we did not complete the corre- 
sponding Fiducial run for N500L34, therefore this run is not 
used for other comparisons in this paper. At the low-mass 
end of Mi, < 10 s Mq, the H2 run has a significantly lower 
number density of galaxies than the Fiducial run. This il- 
lustrates that the H2 model has a greater impact on the 
number density of low-mass galaxies. 

3.3.1 On the overprediction of GSMF 

One of the primary motivations for implementing the H2- 
based SF model was to see if it can remedy the overpredic- 
tion of GSMF at low-mass end due to its natural dependence 
on metallicity as we described in Section 1. In the earlier sec- 
tions, we saw that indeed the liVbased SF model reduces 
the number of low-mass galaxies. However, even with the 
new H2 model, we are still over-predicting the number of 
low-mass objects at M+ = 10 7 ' 8 — 1O 8 ' 6 M0 compared to the 
observational estimate of Gonzalez et al. (2010) at 2 = 6 
(Panel [d]). Therefore the H2 model alone does not seem to 
be able to solve this generic problem of CDM model. Our 
simulations also seem to under-predict the number of mas- 
sive galaxies with > 1O 9 ' 5 M when compared to the 
Gonzalez et al. (2010) observational data at 2 = 6. Jaacks 
et al. (2012a) argued that this difference likely originates 
from the different slope in the M*— SFR relation, where the 
observational estimate was derived by using a crude rela- 
tion from 2 ~ 4 and applied to 2 = 6 assuming that it is 
unchanged. In our simulations, the M*— SFR relation has a 
different slope, and this results in a different slope in the 
GSMF. 

Figure 13d also contains the results of applying the duty 
cycle (DC) corrections (Jaacks et al. 2012b) to our compos- 
ite GSMF both with (dot-dashed line) and without (dotted 
line) accounting for dust extinction. Jaacks et al. (2012b) 
defined the DC as the fraction of time that a galaxy exceeds 
the current HST magnitude limit within a certain A2, and 
characterized it with a sigmoid function as a function of M*. 
According to their result, DC for 2 = 6 makes a relatively 
sharp transition from nearly zero at M* < 1O 7 M0, crosses 
0.5 at M+ ~ 1O 8 M , and to almost unity at M* > 10 9 M Q . 
Using this relation, we can apply a correction for the ob- 
servability of low-mass galaxies, and see the impact of SF 
duty cycle on the observed GSMF. Similarly to the results 
of Jaacks et al. (2012b), our GSMF becomes closer to the 
observational estimate after the DC correction. 



3.3.2 GSMF at 2 = 3 and 

Figure 14 shows the GSMF at 2 = 3 (Panel a) & 2 = 
(Panel b) . Panel (a) is composed of data from the N400L34 
& N600L100 runs, and Panel (b) of N600L100 data. Dashed 
lines represent the Fiducial run, while solid lines represent 
the H2 run. The shaded regions at 2 = 3 represent ob- 
servational estimates of the GSMF at 3 < 2 < 4 (yellow) 
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Figure 13. Panels (a-c): Galaxy stellar mass function for three different box sizes at z = 6, plotted against their respective Fiducial 
run. Panel (d): Composite GSMF of the H2 runs, compared with the Fiducial composite GSMF (black dashed line). Additionally we 
show DC corrections (Jaacks et al. 2012b) with (dot-dashed) and without (dotted) dust correction. The yellow shaded region is the 
observational estimate from Gonzalez et al. (2010). 



and 2 < z < 3 (cyan) from Marchesini et al. (2009), fol- 
lowing Choi & Nagamine (2010). Both sets of simulations 
are in agreement with each other and with observations at 
Ai* > 1O 1O M , which corresponds to Mdm > 10 n ' 5 M Q 
(Figure 8b). A substantial difference between the two SF 
models is again seen in galaxies with M* < 10 10 Mq, but 
this is below the current observable flux limit. 

We may try to understand the discrepancies in the 
GSMF in relation to the SHMR. The difficulty is that the 
SHMR is not per unit volume, hence there is no obvious 
direct link between SHMR and GSMF. Suppose M* in low- 
mass halos is increased uniformly, then the normalization 
of SHMR shifts upwards. At the same time, those galax- 
ies would move from the low-mass bin to higher mass bins, 
and the GSMF would be weighted more towards higher 
mass side. For example, Figure 6b suggests that we are 
producing roughly correct amount of stars in halos with 
Mdm < 10 12 M Q at z = 3, and the agreement in the GSMF 
is not so bad either as shown in Figure 14a. Such a com- 



parison provides a consistency check between SHMR and 
GSMF. 

The shaded region at z — (Panel [b]) is the obser- 
vational estimate from Cole et al. (2001). Our simulations 
agree well with the observation near the knee of GSMF 
(M t ~ 10 10 - 8 - IO^^Mq), but over-predicts at both low 
and high mass end. This over-estimation at A/* > lO n ' 3 M0 
is reflected in the overestimation of the SHMR at A/ to t,200 ~ 
10 13 M Q (Figure 6d), which could be due to a lack of AGN 
feedback in our current simulations. At low-mass end (M* < 
1O 1O ' 5 M ), both models over-predict the GSMF, but the H 2 
run to a lesser degree. 

It is clear that simultaneously matching the SHMR and 
GSMF is not an easy task. We expect the inclusion of AGN 
feedback will assist in improving the high-mass end of our 
simulations at low redshift. The new H2-based SF model 
seems to have improved the relations in regards to the low- 
mass end, but does not fully reconcile the differences. Fur- 
ther improvements to our SN feedback prescriptions (e.g., 
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Figure 14. Panel (a): composite GSMF for z = 3 for both the 
H2 (solid lines) and Fiducial (dashed lines) runs. The data is from 
the N400L34 & N600L100 runs. The shaded regions represent ob- 
servational estimates at 3 < z < 4 (yellow) and 2 < 2 < 3 (cyan) 
from Marchesini et al. (2009). Panel (b): GSMF at z = from 
the N600L100 run. The yellow shaded region is the observation 
from Cole et al. (2001). 

momentum feedback by winds) may be required to achieve 
better agreement with observations. 

3.4 Kennicutt-Schmidt (KS) relationship 



Ideally we would like to reproduce the empirical 
Kennicutt-Schmidt relationship naturally in simulations; 
previously the KS relation was imposed in our SF prescrip- 
tion (Choi & Nagamine 2010) and in many others', therefore 
the results matched the observation well by construction. 
The new H2-based SF model provides two main benefits: 
it is not 'tweaked' to match the KS relation, and it is more 
physically realized in that stars are formed out of cold molec- 
ular gas on a depletion time-scale which is equal to about 
1% of the free-fall time (i.e., with 1% efficiency per free-fall 
time). 

To examine the KS relation, we calculate the column 
density of Hi, H2, and SFRs along the z-axis of each halo in 
our simulation on a uniform grid with a cell size of e 2 . A de- 
tailed description of this process can be found in Nagamine 
et al. (2004a). In Section 2.2, we stated that the H2 model 
was accurate for Z ^ 1O _2 Z0, yet we set our metal floor 
below that at Z — 1O _3 Z0. The model fails at low metallic- 
ities by over-predicting the amount of H2 mass. This is due 
to time-dependent effects being neglected within the ana- 
lytical KMT model (Krumholz & Gnedin 2011). However, 
the over-predicted value may be an accurate estimate of how 
much cold material is present to form stars (Krumholz 2012). 
Therefore we simply assume that the /h 2 value calculated 
by Equation (6) for any gas particle with Z < 1O _2 Z0 is 
actually representative of the amount of cold Hi gas, which 
is available for star formation. 

Figure 15(a-d) shows the KS relation for the N600L100 
simulation at 2 — 6, 3, 1, & 0. Each point in this figure repre- 



sents one pixel on the projected x-y plane, and the contour 
is used to represent all the columns from all halos in the 
simulation box. For each redshift, the panel is broken down 
into three sub-panels: the first being the KS relation for Hi 
gas only, second is for H2 gas only, and lastly for Hi +H2. We 
will refer to these panels as KS-HI, KS-H2, and KS-HIH2, 
respectively. Each panel includes the KS relation given by 
Equation (3) as a solid red line, with the dashed lines rep- 
resenting the range of slope An = ±0.15. 

In KS-HI panel at z = 0, we also overplot the obser- 
vational data from seven nearby spiral galaxies as a blue 
hatched region (Bigiel et al. 2008, hereafter B08). In KS-H2 
panel (2 = 0), we overplot the low surface density observa- 
tions from the Small Magellanic Cloud (SMC) as red trian- 
gles (Bolatto et al. 2011, hereafter known as Bll). Lastly 
in KS-HIH2 panel (2 — 0), we again plot B08 data as a 
blue hatched region, and Bll data from the SMC as a red 
hatched region. 

There are two major processes driving the evolution of 
these plots. The first is gas depletion: as time passes the cold 
molecular gas used to form stars is depleted, and become less 
available at late times. This is most obvious in the decrease 
of E* between 2 — 3, 1, & 0, corresponding to the downturn 
of the SFRD at z < 2 in Figure 7. The second is metal 
enrichment: the longer a simulation runs, the more enriched 
the gas becomes via SF (Figure 10). This process expands 
the distribution of points to the left-hand-side of the plot, 
because higher metal content allows stars to form at lower 
surface densities, as shown in Figure 1. The distribution of 
points broadens from z = 6 to z = 0, indicating greater 
range of metallicities present in the simulations. 

The KS-HIH2 panels include theoretical results from 
the KMT model (Krumholz et al. 2009) to show the same ef- 
fect. The column densities calculated for each pixel represent 
the smoothed value on a relatively large projected scale of e 2 ; 
if we use this value, the model will underpredict /h 2 , since 
it does not account for clumping of the gas on scales below 
our simulation's spatial resolution limit of e = 4.30 /i -1 kpc, 
as well as the path-length along the column. To account for 
this effect, the KMT model multiplies the calculated gas col- 
umn density by a clumping factor "c" (Ehi+h 2 = c x E ca i c ), 
which increases the surface densities to be compared with 
observations. In order to compare the KMT model result 
with our simulation, the theory lines are shifted to lower 
'computed' surface densities (i.e., E ca i c = Ehi+h 2 /c), which 
brings a good agreement between the KMT model results 
and our simulations. In Figure 15, we adopted c = 5. 

For the KS-HI panel at z — 0, we find disagreement 
between simulation and the B08 data (blue hatched region). 
This is a metallicity effect; our simulations do not contain 
enough high-metallicity columns, and the low metallicity 
columns will form stars at higher surface densities in the 
KMT model. 

In the KS-H2 panel at z = 0, our simulation is in 
good agreement with the observations. The E* starts high 
at z — 6, and eases its way to the lower left due to the two 
processes described above. By 2 = 0, the observations lie 
in the center of our simulation data showing a very good 
agreement even for low surface densities. It should be noted 
that we are directly measuring the amount of H2 in our sim- 
ulation, whereas the observers infer this value from the CO 
luminosity. 
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Figure 15. The Kennicutt-Schmidt relation for the N600L100 simulation at z = 6,3, 1, & 0. Each redshift is broken into three panels: 
SFR surface density as a function of Hl(left panel, hereafter KS-HI), H2 (middle panel, KS-H2), and Hi +H2 surface density (right panel, 
KS-HIH2). In each panel, the solid red line represents the empirical KS relation given by Equation (3), and the dashed red lines represent 
the range of slope An = ±0.15. Blue hatched regions in the KS-HI(d) and KS-HIH2(d) panels, and the blue solid contour in KS-H2(d) 
panel are the observations from Bigiel et al. (2008). The red triangles in the KS-H2(d) panel along with the red hatched region in the 
KS-HIH2(d) panel are low-metallicity SMC observations from Bolatto et al. (2011). At z = 3 we have observational data from Rafclski 
et al. (2011) plotted as green circles and black squares. The green circles represent upper limits derived for DLAs, while the black squares 
represent outskirts of LBGs. Lastly in KS-HIH2 panels, black lines represent theoretical results from the KMT model (Krumholz et al. 
2009). In the z = 6 KS-HIH2 panel, the four theory lines correspond to the metallicities log(oZ/Z Q ) = 0.11,-0.69,-1.49,-2.29 from 
left to right, respectively. For z = 3, we have log(cZ/Z ) = 0.47, -0.45, -1.37, -2.29. For 2 = 1, log(cZ/Z ) = 0.77, -0.25, -1.27, -2.29, 
and for z = log(cZ/Zo) = 0.80, —0.23, —1.26, —2.29. Any discreteness of the dotted points at the contour edge is an artifact from our 
plotting procedure. 



In the KS-HIH2 panel at z = 0, we again find a disagree- 
ment between simulation and the B08 data (blue hatched 
region); the bulk of our data is found at slightly higher sur- 
face densities compared to these observations. The data in 
these panels is dominated by Hi, resulting in similar trends 
to the KS-HI panel. 

In the KS-HIH2 panel at z = 3 (Panel [b]), we also over- 
lay the upper limits from damped Lyman-alpha absorbers 
(DLAs) as green circles and outskirts of Lyman-break galax- 
ies (LBGs) as black squares (Wolfe & Chen 2006; Rafel- 
ski et al. 2011). LBGs are considered to be star-forming 
galaxies with moderate median mass of ~ 1O 1O M0, 
therefore are expected to have been enriched to some level. 
Rafelski et al. (2011) find the LBGs in their sample have 
Z » 0.07 - 0.26 Z . Figure 16 shows the KS plot for 
only star-forming columns in our N600L100 simulation with 
Z = 0.07 — 0.26 Zq at z = 3. The observed data points are 
close to the edge of the simulation contour, but there are 
many columns that agree with the observational data. Note 
that it is certainly easier to observe the SFR closer to the 



upper edge of the contour rather than the bottom side of it 
due to the surface brightness dimming. 

Figure 17 further illustrates the metallicity effect by 
separating the KS-HIH2 panel from Figure 15d into three 
different metallicity ranges for the N600L100 run at z = 0. 
We find that the columns with the lowest metallicity (Panel 
[a]) are forming stars at the highest gas surface densities 
for a given E* as expected. Panel (b) brackets Z — 0.2 Zq, 
which is roughly equal to the metallicity of the SMC (Bo- 
latto et al. 2011). Columns in our simulation in this metal- 
licity range agree very well with the observed data (red con- 
tour). In Panel (c) we show columns of higher metallicity 
Z > 0.3 Zq, which is similar to the range of B08 sample 
(0.41 - 0.69 Zq) (Walter et al. 2008). As discussed previ- 
ously, our data does not agree very well with observations 
in this range. There are some points overlapping with the 
observed data, but the majority lie at higher column densi- 
ties than the observed range. This discrepancy is presumably 
caused by different metallicities: the highest metal column 
at z = in the N600L100 run is Z = 1.26 Z , yet the me- 
dian column metallicity at Z > 0.3 Z is Z = 0.41 Zq. This 
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Figure 16. Enlarged region of Figure 15b KS-HIH2 plot. Here we 
only plot star-forming columns with metallicities consistent with 
observations from Rafelski et al. (2011) (Z = 0.07 - 0.26 Zq). 
Green circles represent upper limits derived for DLAs (Wolfe & 
Chen 2006), and black squares represent outskirts of LBGs (Rafel- 
ski et al. 2011). The observed data points are at the upper edge 
of the simulation data contour, but there are many simulated 
columns that overlap with the observed data. 

suggests that our N600L100 run does not contain enough 
high metallicity columns to match these observations. If the 
simulation had more high metallicity columns, then the SF 
would occur more at lower gas surface density, and there 
would be more points overlapping with the B08 data. 

3.5 Hi & H2 column density distribution functions 

One of the best ways to investigate the distribution of Hi 
gas in the Universe statistically is to examine the Hi col- 
umn density distribution function /(Am) (e.g., Nagamine 
et al. 2004b, a; Wolfe et al. 2005; Zwaan & Prochaska 2006; 
Prochaska & Wolfe 2009; Noterdaeme et al. 2009; Pontzen 
et al. 2010; Altay et al. 2011; Yajima et al. 2011; Rahmati 
et al. 2012; Erkal et al. 2012; Noterdaeme et al. 2012). Us- 
ing the Fiducial runs, Nagamine et al. (2010) found that a 
simple self-shielding model with a threshold density (nj^ v = 
6 x 10 -3 cm -3 ) for UVB penetration can reproduce the ob- 
served /(Am) quite well at log Am < 21.5 for z — 3. Yajima 
et al. (2011) later showed the validity of nj^ v value using 
full radiative transfer calculations. However, Nagamine et al. 
(2010) also found that the Fiducial run over-predicts /(Am) 
at log Am > 22, and argued that this might be due to the 
neglect of H2 within the Pressure SF model (Section 2.1.2), 
because then part of Hi would be converted into H2 and 
/(JVhi) would decrease at high Am values. 

Figure 18 compares the column density distributions 
of both Hi and H2 in the H2 and Fiducial runs at z = 
6, 3, 1, & 0. Panels (a) & (b) are composed of N144L10 data, 
while Panels (c) & (d) are composed of N600L100 data. The 
Fiducial run is omitted from panels (c) & (d), because the 
N600L100 Fiducial run did not use the OTUV threshold 
which is necessary to bring the column density distribution 
into agreement with observations at z — 3 (Nagamine et al. 
2010). 



In Panel (a) we see that the H2 run consistently has 
higher amplitude of /(Am) than the Fiducial run due to 
less efficient star formation. At z — 3 (b) however, we find 
that the H2 run has a higher /(Ami) than that of Fiducial 
run at log Ami > 22. This is because the star formation is 
less efficient in the new H2 run, therefore more Hi gas is left 
over in high density regions. In the H2 run the varying SF 
threshold density was higher than the constant nfjf adopted 
in the Fiducial run (Figure 2), and it was also clear from 
Figure 3 that the gas particles are reaching higher densities 
in the H2 run before being heated by SN feedback than in 
the Fiducial run. The f(Nm) results at the lower Am do not 
change between the two runs at this redshift. Panels (c) & 
(d) continue to show the redshift evolution of this relation- 
ship in our simulations. At z = 0, we find that our simula- 
tions over-predict /(Ami) at log Am > 21, over-predict the 
fiHz) at logAm 2 < 21, and under-predict at logAm 2 > 22. 

Therefore the current simulations suggest that it is dif- 
ficult to explain the sharp turn-down of observed /(Am) 
at log Am ~ 22 by the atomic-to-molecular transition, in 
agreement with the conclusions of Erkal et al. (2012). Ad- 
ditionally, Erkal et al. (2012) showed that their simulations 
could be brought into agreement with observations if a re- 
gion of 3 kpc radius around the center of all galaxies was 
removed. This could be another opportunity for AGN feed- 
back to play an important role: if feedback from super mas- 
sive black holes can prevent the formation of high columns, 
then our simulations may come into better agreement with 
observations at A > 10 22 cm -2 . Obviously more refinement 
of feedback models are needed to bring the simulations into 
agreement with the observations of /(Am) and /(Ah 2 ). 

3.6 Resolution studies 

The new H2-based SF model has implicit resolution de- 
pendence. With higher resolution, the simulation resolves 
higher (column) densities (Eq. 9), which yield lower s val- 
ues (Eq. 11) for a given metallicity. Figure 1 illustrates that 
lower s values lead to higher /h 2 (Eq. 6), which increases 
the SFR (Eq. 14). 

To examine the resolution effect, Figure 19 shows 
the KS relation for the N600L10, N400L10, N400L34, & 
N600L100 runs at z — 6. These panels are ordered by reso- 
lution: Panel (a) shows the highest resolution, and Panel (d) 
shows our lowest resolution simulation. In Panels (a) & (b), 
we can examine the resolution effect on the KS plot when 
keeping the box size constant. In general the gas surface 
densities where SF takes place do not change very much, 
but with higher resolution, the points cover a wider range of 
E*. This is an expected result from a higher resolution sim- 
ulation; the additional resolution allows the gas to collapse 
to higher densities, yielding additional shielding which eases 
the transition to H2. 

While Panels (a) & (b) are both from simulations with 
a box size of comoving 10/i _1 Mpc, Panels (c) & (d) arc 
from simulation boxes of 34/i _1 Mpc and 100/i~ Mpc, re- 
spectively. Increasing the box size of a simulation usually 
comes with a price of decreasing the resolution, and it re- 
sults in more higher mass halos and fewer low-mass halos 
(e.g., Thompson & Nagamine 2012). Note that the simula- 
tions shown in Panels (c) & (d) are of lower resolution than 
those in Panels (a) & (b). When comparing the L10 boxes 
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Figure 17. The KS-HIH2 panel from Figure 15d is separated into three different metallicity ranges at z = using the N600L100 run, 
in order to show the metallicity effect on the KS plot. Panel (a) only shows columns with the lowest metallicities. The metallicity range 
in panel (b) brackets the Bolatto et al. (2011) SMC data. Panel (c) shows the highest metal columns, however, the median metallicity of 
the simulated columns are biased towards the lower end of the bracket Z = 0.3 Zq, which is presumably causing the offset between the 
simulation result and the blue hatched observed data. 



with the L34 & L100 boxes, we are actually examining the 
resolution and box size effects simultaneously. Comparing all 
the panels in Figure 19 suggests that our KS results are not 
significantly affected by these resolution effects. The only 
visible effect we see in the figure is that the lower resolution 
results in a thinner contour distribution 

3.6.1 Probability Distribution Function (PDF) 0/H2 
density 

Physical number densities of observed molecular clouds are 
on the order of a few hundred cm -3 , in rough agreement with 
the highest densities achieved in our current cosmological 
simulations. Figure 20 shows the mass-weighted PDF of H2 
number density at the highest densities in our simulations at 
z — 6. As expected, we can see that the peak of the highest 
density region shifts to higher densities as the resolution 
increases; the lowest resolution production run (N600L100) 
has a peak at nn 2 ~ 10 2 cm~ 3 , and our highest resolution 
production run (N400L10) has a peak at nn 2 ~ 10 3 ' 6 cm~ 3 . 
However, the N400L10 run has a slightly different shape 
from the other runs, and the higher resolution N600L10 run 
has a peak at a slightly lower value of riH 2 ~ 10 2,8 cm -3 . 
The exact reason for this different PDF shape is unclear, 
but presumably it was affected by some SF events. 

Earlier, Jaacks et al. (2012a) showed that the Fiducial 
runs do not satisfy the Bate & Burkert (1997) mass resolu- 
tion criteria, even though gas particles in our N400L10 have 
particle masses lower than the typical Jeans mass at z = 6 
by a factor of ~ 1 — 100. This prevents us from explicitly 
resolving the collapse of star forming molecular clouds di- 
rectly, and it is one of the primary reasons for employing a 



sub-grid model for SF using the KMT model. Given that the 
highest densities reached in our simulations is approximately 
equivalent to those of observed giant molecular clouds, we 
consider that the KMT is suitable to use as a sub-grid model 
in our simulation to estimate the H2 mass for star forma- 
tion. In fact, the KMT model is well suited to predict the 
galactic-scale trends in atomic and molecular content rather 
than the structure of individual photo-dissociation regions 
(Krumholz et al. 2008; Kuhlen et al. 2012). 



4 SUMMARY 

We have implemented a new H2-based SF model in our cos- 
mological SPH code GADGET-3 . Previous SF models did 
not consider the formation of H2 , and imposed the KS rela- 
tion in their SF prescriptions. The analytic KMT model has 
provided a computationally inexpensive way to estimate the 
H2 mass fraction in cosmological hydrodynamic simulations, 
which allows us to modify our SF prescription to compute 
the SFR based on H2 density rather than total gas density. 
The model brings a natural dependence of star formation on 
metallicity (in addition to the previous dependence through 
metal line cooling). 

We performed a series of cosmological simulations with 
different box sizes and resolutions, and examined how this 
new H2-based SF model affected the results such as stellar- 
to-halo mass ratio, cosmic SFRD, galaxy specific star forma- 
tion rates, galaxy stellar mass function, Kennicutt-Schmidt 
relationship, and Hi & H2 column density distributions. We 
find that this new H2-based SF model provides many advan- 
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Figure 18. Column density distribution functions of Hi and H2 
at z = 6, 3, 1, & for the H2 and Fiducial runs. Redshifts z = 6 
Sz 3 are from the N144L10 runs, while 2 = 1 & are from the 
N600L100 H 2 runs (N600L100 Fiducial was omitted because it 
lacks the OTUV threshold). The observational data points at 
2 = 3 are from Peroux et al. (2005, black squares), O'Meara 
et al. (2007, magenta squares), Prochaska & Wolfe (2009, green 
triangles), and Noterdaeme et al. (2012, blue circles). Panel (d) 
shows observations from Zwaan & Prochaska (2006), where black 
circles represent the Hi and black triangles represent the H2 col- 
umn density distribution functions. 



tages over previous models, and we summarize our primary 
conclusions below. 

• In the new Eb-based model, each gas particle has dif- 
ferent SF threshold densities based on its metallicity (Fig- 
ure 2). We have shown that the new SF threshold densities 
(i.e., metallicity-dependent density required for H2 forma- 
tion) are higher than the constant threshold density used in 
the Fiducial run, which results in overall decrease of SFRD 
(Figure 7) in the new model. Decrease of star formation leads 
to weaker feedback effects subsequently. The need for suffi- 
cient shielding from radiation field for H2 formation results 
in lower SFR, causing a gas reservoir to build up. Conse- 
quently, SF starts later than in the Fiducial run, and the 
peak of SFRD has slightly shifted to a lower redshift. But 
both runs are still compatible with the observed range of 
SFRD in the Lilly-Madau diagram. 

• The H2 run is able to successfully reproduce the SHMR 
at z = 3 & 6 for lower mass halos with M to t,200 < 1O 12 M0 
(Figure 6). The Fiducial run with previous SF model signifi- 
cantly overpredicts SHMR at the same mass range, therefore 
the H2 run provides a significant improvement on this as- 
pect. Since the SN feedback model was kept the same in the 
two runs, this improvement was purely driven by the change 
in the SF model, rather than the feedback. 



back in our current simulations. This is connected with the 
over-prediction of GSMF at the high-mass end in our simu- 
lations. 

• The sSFRs of galaxies in the H2 and Fiducial runs are in 
rough agreement with observations (Figure 11), and they de- 
crease systematically with decreasing redshift. At z — 6, the 
H 2 run have higher sSFR for galaxies with M* < 10 10 M Q , 
but this is due to the fact that the galaxies with same M* 
reside in higher mass halos in the H2 runs than in the Fidu- 
cial run (Figure 8). At later times, this difference becomes 
much smaller and the two models are in rough agreement 
with one another. 

However, the median sSFR of simulated galaxies with 
M* = 1O 1O M0 does not behave consistently with the ob- 
servations as a function of redshift (Figure 12). Both our 
Fiducial and H2 runs predict gradually declining sSFR 
with decreasing redshift, similarly to other previous mod- 
els (Bouche et al. 2010; Dave et al. 2011; Weinmann et al. 
2011). Krumholz & Dekel (2011) argued that taking the 
metallicity-dependence of H2 formation would help to recon- 
cile the discrepancy, however, even with our new H2-based 
SF model, our simulations do not produce the plateau of 
sSFR at z > 2. The general agreement between multiple 
different simulations and semi-analytic models of galaxy for- 
mation suggest that the ACDM model predicts a general 
decline in the sSFR of galaxies of a given mass, contrary to 
observations. However, we note that none of these simula- 
tions included the effect of AGN feedback. 

• We find that the H2-based SF model produces signifi- 
cantly fewer galaxies at M* < 1O 8 M0 compared to the Fidu- 
cial run at 2 = 6 (Figure 13d). Even after this reduction, the 
faint-end slope of GSMF in the H2 run is still steeper than 
what has been observationally estimated at z = 6. Employ- 
ing duty cycle corrections following Jaacks et al. (2012b) 
brings the GSMF closer to observations. 

At z = 3 we find that our simulations are in good agree- 
ment with observed GSMFs at M* > 10 10 M Q , consistent 
with our previous finding in Choi & Nagamine (2010). At 
the lower masses of M* < 1O 1O M0, again the H2 model 
produces fewer number of low-mass galaxies relative to the 
Fiducial run. At the moment, the flux limit of GSMF data 
is M* ~ 1O 1O M even with the deepest HST imaging, 
and there is no good data below this limit. Galaxies with 
M* < 1O 1O M correspond to halos with Mdm < 1O 12 M , 
and in this regime the new H2 run agrees with the observa- 
tional estimate of SHMR much better than the Fiducial run. 
For this reason, we expect that the H2 run would match the 
observations of GSMF better in the future at M* < 10 10 M Q . 

Finally at z — 0, we find that our simulations over-predict 
the GSMF at both low and high-mass end. The deviation at 
the low-mass end seems smaller than at the high-mass end, 
however, since this is a log-log plot, the actual deviation is 
greater at the low-mass end. Further improvement in our 
feedback prescriptions (e.g., momentum-driven feedback by 
SN and AGN) may be needed to reconcile these differences. 

• We find that the new H2-based SF model naturally pro- 
duces the empirical Kennicutt-Schmidt relationship without 
the need for 'tweaking' the parameters of the SF model. The 
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Figure 19. Similarly to Figure 15, we plot the KS relation for the N600L10, N400L10, N400L34, & N600L100 simulations at z = 6 to 
examine the resolution and box size effects. See text for detailed discussions. 



most significant discrepancy between the H2 run and obser- 
vation can be seen against the nearby spiral galaxy data of 
Bigiel et al. (2008). It seems that the H2 run does not con- 
tain enough high-metallicity columns to match observations 
of nearby spiral galaxies, and the same trend can be seen in 
the galaxy mass-metallicity plot (Figure 10). These discrep- 
ancies indicate that our current simulations might have too 
much low-metallicity gas in massive galaxies at z < 2, which 
is also indicated by Figure 9d. However the H2 run is able to 
match the observations of DLAs and LBGs at z = 3, as well 
as the observations of the low-metallicity SMC by z — 0. 

• As for the hydrogen column density distribution func- 
tion, we find that the new H2 model did not improve the 
agreement with observation at log Nh > 21.6 at 2 = 3, and 
we still over-predict f(Nni) similarly to the previous simula- 
tions. Erkal et al. (2012) also concluded that the atomic-to- 
molecular transition alone could not account for the down- 
turn in /(-/Vhi) at log Nh > 21. At z — 0, our simulations do 
not agree with the observational data of Zwaan & Prochaska 
(2006), and further refinement of star formation and feed- 
back models are needed. 

As for our future plan, we intend to improve our simu- 
lations on a few fronts, given the problems that we observed 
in this paper. Since the /h 2 calculated in the KMT model 
depends on gas metallicity, we need to account for the metal 
diffusion in the ISM more accurately (e.g. Shen et al. 2010). 
Our current SPH code does not allow for particles to share 
their metal content with one another, and we plan to imple- 
ment and explore the effects of metal diffusion in the near 
future. Finally, as a comparison to the H2-based SF model, 
we also plan to develop a turbulence-based SF model and 



and explore the differences between the two approaches to 
star formation. 
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